{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5678b63a-3352-4256-8691-d9578a292fa9",
   "metadata": {},
   "source": [
    "## Estimate structural parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "56a0e031-ed12-4bb1-9820-7183d7a4c0e4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "from fn_simulate import *\n",
    "\n",
    "np.set_printoptions(precision=4)\n",
    "\n",
    "# randomization variables for simulation\n",
    "#import secrets\n",
    "#secrets.randbits(128) \n",
    "seed = 84708628852577542122415832887111152745\n",
    "\n",
    "rng = np.random.default_rng(seed)\n",
    "\n",
    "n_draws = 1000\n",
    "\n",
    "param_random = {'n_draws' : n_draws, \n",
    "                'seed' : seed}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ca330fd5-a3cf-4d09-a8af-b0820508d2f9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# calibrated parameters\n",
    "\n",
    "    # income\n",
    "mu_lny = 9.713\n",
    "sig_lny = 0.497\n",
    "\n",
    "mu_eps = 0\n",
    "\n",
    "    # prices\n",
    "pi_s = 1\n",
    "pi_q = 0\n",
    "tau = 0.075\n",
    "\n",
    "gamma = 0.424\n",
    "\n",
    "param_given = {'mu_lny' : mu_lny, \n",
    "               'sig_lny' : sig_lny,\n",
    "               'mu_eps': mu_eps,\n",
    "               'pi_s' : pi_s, \n",
    "               'pi_q' : pi_q, \n",
    "               'tau' : tau,\n",
    "               'gamma': gamma}\n",
    "\n",
    "\n",
    "# estimated parameters\n",
    "# initial: [ 1.467219  0.110516  0.92719   0.611335 -7.924372  0.299971 -0.139612 -0.010366]\n",
    "# results:  [ 1.466967  0.110715  0.927835  0.612176 -7.919992  0.299912 -0.140181 -0.009912]\n",
    "\n",
    "# initial:\n",
    "pi_n = 1.467219\n",
    "pi_nq = 0.110516\n",
    "\n",
    "theta = 0.927190\n",
    "alpha = 0.611335\n",
    "rho = -7.924372\n",
    "\n",
    "sig_eps = 0.299971\n",
    "eps_cutoff = -0.139612\n",
    "\n",
    "cor_lny_eps = -0.010366\n",
    "\n",
    "\n",
    "param_choice = np.array( \n",
    "    [pi_n, pi_nq, \n",
    "     theta, alpha, rho,\n",
    "     sig_eps, eps_cutoff, cor_lny_eps] \n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "65173b86-59ad-4ad6-97e0-35cc4b5b5387",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- 7.176711797714233 seconds ---\n",
      "beta_A_2to3:  -0.15696830315310745\n",
      "beta_B_2to3:  0.4119744870325511\n",
      "PA_2to3:  0.21208907741251326\n",
      "PB_2to3:  0.12089077412513255\n",
      "exp_n_AB:  0.21271942877668723\n",
      "exp_nq_AB:  0.11103227872559318\n",
      "cor_lny_n3:  0.12864374079269372\n",
      "cor_lny_lnq_AB:  0.29604023673643337\n"
     ]
    }
   ],
   "source": [
    "# trial\n",
    "start_time = time.time()\n",
    "moments_simu = simulate_moments(param_choice, param_given, param_random, disp = False)\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "\n",
    "print('beta_A_2to3: ', moments_simu[0])\n",
    "print('beta_B_2to3: ', moments_simu[1])\n",
    "print('PA_2to3: ', moments_simu[2])\n",
    "print('PB_2to3: ', moments_simu[3])\n",
    "print('exp_n_AB: ', moments_simu[4])\n",
    "print('exp_nq_AB: ', moments_simu[5])\n",
    "print('cor_lny_n3: ', moments_simu[6])\n",
    "print('cor_lny_lnq_AB: ', moments_simu[7])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bccf8ec9-4482-4e72-8528-7f09b1a5f5e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# moments in the real data\n",
    "\n",
    "# moments observed in real data\n",
    "beta_A_2to3_real = -0.157  \n",
    "beta_B_2to3_real = 0.412 \n",
    "\n",
    "PA_2to3_real = 0.212\n",
    "PB_2to3_real = 0.121\n",
    "\n",
    "exp_n_AB_real = 0.214\n",
    "exp_nq_AB_real = 0.111\n",
    "\n",
    "cor_lny_n3_real = 0.128\n",
    "cor_lny_lnq_AB_real = 0.296\n",
    "\n",
    "moments_real = np.array( \n",
    "    [beta_A_2to3_real, beta_B_2to3_real, PA_2to3_real, PB_2to3_real, \n",
    "     exp_n_AB_real, exp_nq_AB_real, cor_lny_n3_real, cor_lny_lnq_AB_real] \n",
    ")\n",
    "\n",
    "# weights: 1/moments_real^2 as in Baudin, de la Croix, & Gobbi (2015, AER)\n",
    "W = np.diag( 1 / np.square(moments_real) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8c43c491-3281-4709-ac76-a44f940e39f7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Display for each simulation: param_choice and value of the objective function\n",
    "disp = False\n",
    "\n",
    "# Define the objective function to be minimized\n",
    "def objective_function(param_choice): # choice variables should be a numpy array\n",
    "    \n",
    "    return simulate_objective(simulate_moments, param_choice, param_given, param_random, moments_real, W, disp = disp)\n",
    "\n",
    "# Initial guess for the choice parameters\n",
    "param_choice0 = param_choice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "74229357-72f3-4862-9d46-3c959a4d73ef",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.000057\n",
      "         Iterations: 199\n",
      "         Function evaluations: 359\n",
      "--- 2567.7739691734314 seconds ---\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(precision=4)\n",
    "\n",
    "#import time\n",
    "start_time = time.time()\n",
    "\n",
    "# Perform the optimization\n",
    "result = minimize(objective_function, param_choice0, method = 'Nelder-Mead', \n",
    "                  options = {'disp':True, 'fatol':0.0001, 'maxiter':1000}\n",
    "                 )\n",
    "\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "2cb4f3df-2b35-4dba-b5bb-ebde267b6109",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the results\n",
    "param_choice_opt = result.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bda0fe6a-0b96-42f6-9e5a-289a2b5b6aca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[beta_A_2to3, beta_B_2to3, PA_2to3, PB_2to3,exp_n_AB, exp_nq_AB, cor_lny_n3, cor_lny_lnq_AB]\n",
      "[ 1.466967  0.110715  0.927835  0.612176 -7.919992  0.299912 -0.140181\n",
      " -0.009912]\n"
     ]
    }
   ],
   "source": [
    "np.set_printoptions(precision=6)\n",
    "print('[beta_A_2to3, beta_B_2to3, PA_2to3, PB_2to3,exp_n_AB, exp_nq_AB, cor_lny_n3, cor_lny_lnq_AB]')\n",
    "print(param_choice_opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d364941b-24f7-4468-b9d2-c0128dc9d450",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_A_2to3:  -0.15686103860272904\n",
      "beta_B_2to3:  0.4119390424079334\n",
      "PA_2to3:  0.21208907741251326\n",
      "PB_2to3:  0.12089077412513255\n",
      "exp_n_AB:  0.21266544572412097\n",
      "exp_nq_AB:  0.11121397096328621\n",
      "cor_lny_n3:  0.12845500745776375\n",
      "cor_lny_lnq_AB:  0.2960323331376951\n"
     ]
    }
   ],
   "source": [
    "# simulate using the estimated parameter values\n",
    "\n",
    "pi_n = 1.466967\n",
    "pi_nq = 0.110715\n",
    "\n",
    "theta = 0.927835\n",
    "alpha = 0.612176\n",
    "rho = -7.919992\n",
    "\n",
    "sig_eps = 0.299912\n",
    "eps_cutoff = -0.140181\n",
    "\n",
    "cor_lny_eps = -0.009912\n",
    "\n",
    "param_choice = np.array( \n",
    "    [pi_n, pi_nq, \n",
    "     theta, alpha, rho,\n",
    "     sig_eps, eps_cutoff, cor_lny_eps] \n",
    ")\n",
    "\n",
    "moments_simu = simulate_moments(param_choice, param_given, param_random, disp = False)\n",
    "\n",
    "print('beta_A_2to3: ', moments_simu[0])\n",
    "print('beta_B_2to3: ', moments_simu[1])\n",
    "print('PA_2to3: ', moments_simu[2])\n",
    "print('PB_2to3: ', moments_simu[3])\n",
    "print('exp_n_AB: ', moments_simu[4])\n",
    "print('exp_nq_AB: ', moments_simu[5])\n",
    "print('cor_lny_n3: ', moments_simu[6])\n",
    "print('cor_lny_lnq_AB: ', moments_simu[7])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:base] *",
   "language": "python",
   "name": "conda-base-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
